Library loading and Set up

First, Set up the random seed and sanitize the working environment.

rm(list = ls())
set.seed(2379)

Then, load the libraries.

library(dplyr)
library(ggplot2)
library(gridExtra)
library(tidyverse)
library(ROCR)
library(knitr)
library(tibble)
library(plotly)
library(lime)
library(rpart)
library(rpart.plot)
library(dendextend)
library(vtreat)
library(class)
library(grDevices)
library(factoextra)
library(fpc)

Load Data

Assign the dataset to the variable ‘YT’ after reading the csv file and importing it into R.

YT <- read.csv("./youtube_UTF_8.csv", header = TRUE, encoding = "UTF-8")

Data Cleaning

# Replace all the missing values with NA
YT[YT == ""] <- NaN

YT <- YT %>%
  mutate(
    allunemployment = (Unemployment.rate / 100) * Population,
    unemployment_urban = (Unemployment.rate / 100) * Urban_population,
    educationall = (Gross.tertiary.education.enrollment.... / 100) * Population,
    educationurban = (Gross.tertiary.education.enrollment.... / 100) * Urban_population
  )

YT <- subset(YT, select = -c(Youtuber, Longitude, Latitude, Title, country_rank, Abbreviation, channel_type_rank,Gross.tertiary.education.enrollment....,Unemployment.rate,video_views_rank))


YT$created_year <- as.character(YT$created_year)
YT$created_date <- as.character(YT$created_date)


# replace all the "nan" values in categorical variables to "missing"
YT[] <- lapply(YT, function(x) {
  if (is.character(x)) {
    x[x == "nan"] <- "missing"
  }
  return(x)
})

# delete the rows that have missing values in all numerical variables except Longitude and Latitude
YT <- YT[!(rowSums(is.na(YT[, -c(1, 2)])) == ncol(YT) - 2), ]

# Change missing values with the median value in the certain columns
median_re <- function(YT) {
  nucol <- setdiff(names(YT), c("created_year", "created_date", "rank"))
  
  for (col in nucol) {
    if (is.numeric(YT[[col]])) {
      YT[[col]][is.na(YT[[col]])] <- median(YT[[col]], na.rm = TRUE)
    }
  }
  
  return(YT)
}
YT <-  median_re(YT)

cand.columns <- names(YT)[names(YT) != 'subscribers']
str(YT)

## 'data.frame':    995 obs. of  22 variables:
##  $ rank                            : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ subscribers                     : int  245000000 170000000 166000000 162000000 159000000 119000000 112000000 111000000 106000000 98900000 ...
##  $ video.views                     : num  2.28e+11 0.00 2.84e+10 1.64e+11 1.48e+11 ...
##  $ category                        : chr  "Music" "Film & Animation" "Entertainment" "Education" ...
##  $ uploads                         : int  20082 1 741 966 116536 0 1111 4716 493 574 ...
##  $ Country                         : chr  "India" "United States" "United States" "United States" ...
##  $ channel_type                    : chr  "Music" "Games" "Entertainment" "Education" ...
##  $ video_views_for_the_last_30_days: num  2.26e+09 1.20e+01 1.35e+09 1.98e+09 1.82e+09 ...
##  $ lowest_monthly_earnings         : num  564600 0 337000 493800 455900 ...
##  $ highest_monthly_earnings        : num  9.0e+06 5.0e-02 5.4e+06 7.9e+06 7.3e+06 ...
##  $ lowest_yearly_earnings          : num  6.8e+06 4.0e-02 4.0e+06 5.9e+06 5.5e+06 ...
##  $ highest_yearly_earnings         : num  1.08e+08 5.80e-01 6.47e+07 9.48e+07 8.75e+07 ...
##  $ subscribers_for_last_30_days    : num  2e+06 2e+05 8e+06 1e+06 1e+06 2e+05 2e+05 2e+05 1e+05 6e+05 ...
##  $ created_year                    : chr  "2006" "2006" "2012" "2006" ...
##  $ created_month                   : chr  "Mar" "Mar" "Feb" "Sep" ...
##  $ created_date                    : chr  "13" "5" "20" "1" ...
##  $ Population                      : num  1.37e+09 3.28e+08 3.28e+08 3.28e+08 1.37e+09 ...
##  $ Urban_population                : num  4.71e+08 2.71e+08 2.71e+08 2.71e+08 4.71e+08 ...
##  $ allunemployment                 : num  73239992 48251210 48251210 48251210 73239992 ...
##  $ unemployment_urban              : num  25247290 39787465 39787465 39787465 25247290 ...
##  $ educationall                    : num  3.84e+08 2.90e+08 2.90e+08 2.90e+08 3.84e+08 ...
##  $ educationurban                  : num  1.32e+08 2.39e+08 2.39e+08 2.39e+08 1.32e+08 ...

Firstly, the missing values are replaced with ‘NaN’ and irrelevant columns are removed. ‘Missing’ is the default value for categorical variables containing the character ‘nan’. In order to identify interactions such as unemployment rates multiplied by population metrics and tertiary education enrollment rates multiplied by population, feature engineering is used to construct new variables based on current columns. Particularly, rows with entirely missing numerical data, apart from Longitude and Latitude, are eliminated. The data integrity is then maintained by using medians to impute missing numerical values. Data cleaning is vital for maintaining the dataset’s quality and dependability for analytical purposes and for preparing it for later machine learning modelling.

Target Variable

In the analysis, the target (response) variable selected is a newly indicated column subscribersup, which indicates the number of YouTube channel subscribers. This variable has been carefully crafted to create a distinct difference between channels that have a large and low subscriber base. Based on the median value, the dataset is separated into two groups: the majority of channels have subscribers below 17,700,000, while the top 5% of the data are channels with subscribers at or above this level. Channels with a lower subscriber count are labeled [1, 1.77e+07], while channels with a greater subscriber count are labeled [1.77e+07, Inf]. The target variable is designed to offer practical insights into audience engagement, content strategy, and advertising prospects as a focus point of analysis for enhancing YouTube channel performance and impact.

YT <- YT %>%
  mutate(subscribersup = cut(subscribers, breaks = c(0, 17700000, Inf), include.lowest=T, right=F)) %>% 
  
  filter(!is.na(subscribersup))
table(YT$subscribersup)

## 
##   [0,1.77e+07) [1.77e+07,Inf] 
##            490            505

YT$subscribersupbool <- YT$subscribers >= median(YT$subscribers)
target <- 'subscribersup'
pos.label <- "[1.77e+07,Inf]"
YT <- YT %>%
  select(-subscribers)

By examining the proportions of outcome variables selected, it can be observed that the outcome is distributed equally. Building reliable and unbiased machine learning models requires balanced representation of both classes: low-subscriber and high-subscriber, which is accomplished by an equal distribution. Balanced class is essential in the context of classification tasks since it avoids the model from being biased towards predicting the majority class, ensuring fair evaluation and trustworthy insights into the reasons influencing high subscriber counts.

Feature Variable

In order to improve the effectiveness and accuracy of machine learning models, feature selection is a vital step. To reduce complexity and improve model performance, it involves selecting relevant features from a dataset (Khalid, Khalil, & Nasreen, 2014). The correlation coefficient stands out among the numerous techniques as a useful statistical measure, illuminating the strength and direction of correlations between variables. This coefficient shows to be an essential tool in the field of machine learning, assisting in the discovery and choice of the features that have the greatest impact on the target variable (Hsu & Hsieh, 2010).

Benefits of Correlation Coefficient Feature Selection:

1. Improve Model Performance: The model’s ability to recognize and generalize patterns is strengthened by removing unnecessary or redundant information. As a result, this improvement produces more accurate predictions on unobserved data, which is essential for guaranteeing the model’s applicability and trustworthiness in the actual world (Tang, Alelyani, & Liu, 2014).
2. Enhance Computational Efficiency: Machine learning methods are less computationally complex when they operate on a reduced feature set. This decrease in complexity speeds up not just the training procedure but also the inference process. The ability to process enormous volumes of data quickly is crucial in the fast-paced world of current technology, making this feature selection consideration of crucial significance (Elssied, Ibrahim, & Osman, 2014).
3. Interpretability and Insights: Data scientists can better understand model interpretation and reveal complex correlations between features and the target variable using a smaller feature set. Understanding these relationships not only promotes better model comprehension but also provides data scientists with practical knowledge. The findings of the model are made more accessible and beneficial for stakeholders because to this interpretative clarity, which is crucial in decision-making processes (Karatza, Dalakleidi, Athanasiou, & Nikita, 2021).

The correlation coefficient is employed for feature selection in the following data modelling process due to its proven efficiency. The advantages listed above, such as improved model performance, increased computational efficiency, and improved interpretability, make it clear why this method is essential for use in data science and machine learning. The correlation coefficient is used to ensure that machine learning models are accurate, efficient, and insightful, optimizing their influence across a variety of domains.

# Identify categorical and numerical variables
YTclass <- setdiff(colnames(YT), c(target,'subscribersup'))

YTCV <- YT %>%
  select(all_of(YTclass)) %>%
  select_if(~ is.factor(.) || is.character(.)) %>%
  names()

YTNV <- YT %>%
  select(all_of(YTclass)) %>%
  select_if(~ is.numeric(.)) %>%
  names()

YTNV_cluster <- YT[YTNV]
cat("There are ",length(YTCV)," categorical variables and ",length(YTNV), "numerical variables.")

## There are  6  categorical variables and  15 numerical variables.

The dataset YT’s is being divided into two categories: categorical and numerical variables. Categorical variables are identified based on their factor or character data types whereas numerical variables are recognized based on their numeric data type.

Convert Categorical variables to numerical variables used for creating the single variable model, but still keep the original variable for future multivariate model

for (cV in YTCV) {
  catcount <- table(YT[, cV])
  countcv <- paste("Numerical", cV, sep = '_')
  YT[countcv] <- catcount[match(YT[, cV], names(catcount))]
}

YTnewclass <- setdiff(colnames(YT), target)
YTCV <- YT %>%
  select(all_of(YTnewclass)) %>%
  select_if(~ is.factor(.) || is.character(.)) %>%
  names()

YTNV <- YT %>%
  select(all_of(YTnewclass)) %>%
  select_if(~ is.numeric(.)) %>%
  names()

YTNVdata <- YT[sapply(YT, is.numeric) | sapply(YT, is.integer)]
YTCVdata <- YT[sapply(YT, is.factor) | sapply(YT, is.character)]

The YT dataset is subjected to a transformation method. The approach initially counts the instances of distinct values for each categorical variable in YTCV and generates a new numerical variable to reflect these counts. In order to enable numerical analysis and modelling, the transformation effectively converts categorical variables into numerical variables.

Training and testing model

Using 60/20/20 method to separate the youtube dataset into Training, Calibration, and Testing data.

YTsaperate <- runif(nrow(YT))
YTtrain <- YT[YTsaperate <= 0.6, ]
YTCal <- YT[YTsaperate > 0.6 & YTsaperate <= 0.8, ]
YTtest <- YT[YTsaperate > 0.8, ]

cat("The training set has", nrow(YTtrain), "observations, the calibration set has", nrow(YTCal)," observations,and the test set has", nrow(YTtest), "observations.")

## The training set has 613 observations, the calibration set has 179  observations,and the test set has 203 observations.

The ‘runif’ function is used to do the splitting of the data at random, giving each row an equal chance of being assigned to each sets.

Null model

A Null Model is created, which is a simple approach that makes direct assumptions about the predictive power of the model beyond the positive class’s base rate.

YTNullm <- sum(YT$subscribersupbool) / nrow(YTtrain)

YTAUC <- function(predcol, outcol, pos=pos.label) {
  perf <- performance(prediction(predcol, outcol==pos),'auc')
  as.numeric(perf@y.values)
}

YTlog <- function(outCol, predCol, pos=pos.label) {
  sum(ifelse(outCol==pos, log(predCol), log(1-predCol)))
}

baseRateCheck <- YTlog(YTCal[,target],  sum(YTCal[,target]==pos.label)/length(YTCal[,target]) )

modelsummary <- data.frame(
  Name = 'Null Model',
  Type = 'univariate',
  AUC = YTAUC(rep(YTNullm, nrow(YTtest)), YTtest[, target])
)

kable(modelsummary)

Single Variable

Create a function that will be used for model building:

mkPredC <- function(outCol,varCol,appCol, pos=pos.label) {
  pPos <- sum(outCol==pos)/length(outCol)
  naTab <- table(as.factor(outCol[is.na(varCol)]))
  pPosWna <- (naTab/sum(naTab))[pos]
  vTab <- table(as.factor(outCol),varCol)
  pPosWv <- (vTab[pos,]+1.0e-3*pPos)/(colSums(vTab)+1.0e-3)
  pred <- pPosWv[appCol]
  pred[is.na(appCol)] <- pPosWna
  pred[is.na(pred)] <- pPos
  pred
}

mkPredN <- function(outCol,varCol,appCol) {
  cuts <- unique(as.numeric(quantile(varCol, probs=seq(0, 1, 0.1), na.rm=T)))
  varC <- cut(varCol, cuts)
  appC <- cut(appCol, cuts)
  mkPredC(outCol, varC, appC)
}

AUC scores are computed for each variable on both the training and calibration sets. Variables having AUC values over 0.5, showing predictive power superior to chance, are chosen for further analysis.

The AUC result of single variable model is sorted and showned.

AUCYTNV <- data.frame(Pred = character(), Type = character(), TrainAUC = numeric(), CalAUC = numeric())
for (v in YTNV) {
  predii <- paste('pred', v, sep = '.')
  YTtrain[, predii] <- mkPredN(YTtrain[, target], YTtrain[, v], YTtrain[, v])
  YTtest[, predii] <- mkPredN(YTtrain[, target], YTtrain[, v], YTtest[, v])
  YTCal[, predii] <- mkPredN(YTtrain[, target], YTtrain[, v], YTCal[, v])
  trainsetAUC <- YTAUC(YTtrain[, predii], YTtrain[, target])
  
  if (trainsetAUC >= 0.5) {
    CalsetAUC <- YTAUC(YTCal[, predii], YTCal[, target])
    AUCYTNV <- rbind(AUCYTNV, data.frame(Predictor = predii, Type = 'univariate', TrainAUC = trainsetAUC, CalAUC = CalsetAUC))
  }
}

AUCYTNV <- AUCYTNV %>%
  filter(!grepl("pred of subscribers", Predictor)) %>%
  mutate(Differ = abs(TrainAUC - CalAUC)) %>%
  arrange(Differ) %>%
  select(-Differ) 
kable(AUCYTNV, digits = 4, align = 'lccc')

ggplotly(ggplot(YTCal) +
  geom_density(aes(x=pred.rank, color=as.factor(subscribersup))) +
  labs(color = 'subscribers') +
  theme_minimal())

Based on the double density plot, it can be observed that there is an absence of overlapping between two categories suggests that the pred.rank effectively distinguishes between low subscriber counts channels and high subscriber count channels.

Log likelihood is computed for each single variable in the calibration set in order to evaluate how well the model fits. It is essential that the variables chosen can improve the prediction in order to make the best choices feasible. The variable is statistically significant if the log likelihood is close to the base rate check value.

selPredVars <- character(0)
selVars <- character(0)
minStep <- baseRateCheck

YTlogNV <- data.frame(Pred = character(0), 'Log likelihood' = numeric(0))


for (v in YTNV) {
  predii <- paste('pred', v, sep = '.')
  likehoodlog <- sum(ifelse(YTCal[, target] == pos.label, log(YTCal[, predii]), log(1 - YTCal[, predii])))
  
  if (likehoodlog > minStep) {
    selPredVars <- c(selPredVars, predii)
    selVars <- c(selVars, v)
    YTlogNV <- rbind(YTlogNV, data.frame(Pred = predii, loglike = likehoodlog))
  }
}

kable(YTlogNV, digits = 4, align = 'lc')

select_variable_log <- as.character(YTlogNV$Pred)

cormap <- YTNVdata %>%
  cbind(subscriberupbool = YT$subscribersupbool) %>%
  cor() %>%
  as.table() %>%
  as.data.frame() %>%
  filter(Var1 == "subscriberupbool") %>%
  filter(Var2 != "subscriberupbool") %>%
  filter(Var2 != "subscribers") %>%
  distinct(Var2, .keep_all = TRUE) %>%
  filter(Freq >= 0.05)


select_variable_corelation <- as.character(cormap$Var2)
select_variable_ODDcorelation <- as.character(cormap$Var2)[seq(1, length(as.character(cormap$Var2)), by = 2)]
select_variable_EVENcorelation <- as.character(cormap$Var2)[seq(2, length(as.character(cormap$Var2)), by = 2)]
select_variable_log <- as.character(YTlogNV$Pred)
kable(cormap)

Another technique employed is examining the correlation between the predictor variables and the target variables. A significant relationship is shown by a correlation coefficient that is higher than a certain threshold, 0.05 in this case. Additionally, this technique is further refined by separating the selected variables into two group: odd-indexed and even-indexed variable for more understanding of the relationships between the predictor variables and the target variable. As a result, there are total of eleven variables chosen as feature variables as listed above.

Multivariate models

Multimodellikelog <- data.frame(name = character(), 'model Log likelihood' = numeric(), 'Null log likelihood' = numeric(0))
Multimodelauc <-  data.frame(Pred = character(), Type = character(), TrainAuc = numeric(), CalAUC = numeric())
YTmodeleval <- function(namea, predmodel, trueV, do.print = TRUE, threshold = 0.5, pos = pos.label, selection = "all"){
  evalmodelauc <- YTAUC(predmodel, trueV)
  if (do.print == TRUE) {
    if (selection %in% c("all", "conma")) {
      YTCMtable <- table(pred = predmodel > threshold, true.values = trueV)
      
      TP <- YTCMtable[2, 2]
      FP <- YTCMtable[1, 2]
      FN <- YTCMtable[2, 1]
      
      precision <- TP / (TP + FP)
      recall <- TP / (TP + FN)
      
      cat("Precision:", precision, "\nRecall:", recall, "\nConfusion matrix for", namea, "\n")
      print(kable(YTCMtable))
    }
    
    if (selection %in% c("all", "loglike")){
        Nullpre <- mean(trueV == pos.label)
        YTmultilog <- YTlog(trueV, predmodel)
        nullmultilog <- YTlog(trueV, Nullpre)
        rowlog<- data.frame(name = namea, 'model Log likelihood' = YTmultilog, 'Null log likelihood' = nullmultilog)
        Multimodellikelog <<- rbind(Multimodellikelog, rowlog)
        cat("Model Log Likelihood is", YTmultilog, "and null Log Likelihood is", nullmultilog, ".\n")
    }
    
    if (selection %in% c("all", "roc")) {
      perf <- performance(prediction(predmodel, trueV), "tpr", "fpr")
      YTROCMulti <- data.frame(
        FPR = unlist(perf@x.values),
        TPR = unlist(perf@y.values),
        threshold = unlist(perf@alpha.values)
      )
      
      print(ggplot(YTROCMulti, aes(x = FPR, y = TPR)) +
        geom_line(color = "blue") +
        geom_abline(intercept = 0, slope = 1, linetype = 'dashed', color = "red") +
        labs(x = "FPR", y = "TPR", title = paste("ROC Curve for", namea)) +
        theme_minimal())
    }
    if (selection %in% c("all", "Predpplot")) {
        PredPplot <- data.frame(predi = predmodel, truevalue = trueV)
        
        plot <- ggplot(PredPplot, aes(x = predi, color = truevalue)) +
          geom_density() +
          labs(x = "Predicted probability", title = paste("Prediction distribution plot for", namea)) +
          theme_minimal()
        
        print(plot)
      }

    if (selection == "all" || selection == "Prereplot") {
          YTprec <- performance(prediction(predmodel, trueV), measure="prec")
          YTrec <- performance(prediction(predmodel, trueV), measure="rec")
          precision <- (YTprec@y.values)[[1]]
          YTxprec <- (YTprec@x.values)[[1]] 
          recall <- (YTrec@y.values)[[1]]
          rocFrame <- data.frame(threshold=YTxprec, precision=precision, recall=recall)
          print(ggplot(rocFrame, aes(x=threshold)) +
              geom_line(aes(y=precision, color="precision")) +
              geom_line(aes(y=recall, color="recall")) +
              coord_cartesian(xlim = c(0, 1), ylim = c(0, 1)) +
              labs(x = "threshold", y = "percentage", title = paste("precision and recall plot for", namea)) + 
              theme_minimal())
        }

  }
  evalmodelauc
}   

YTevalmodelauc <- function(namea, type, name, predmodelfortrain, truefortrain, predmodelforcal, trueforcal, plotname){
  mtrainauc <- YTmodeleval(namea, predmodelfortrain, truefortrain, do.print = FALSE, selection = plotname)
  mcalauc <- YTmodeleval(namea, predmodelforcal, trueforcal, do.print = TRUE, selection = plotname)
  mresultauc <<- data.frame(Pred = namea, Type = type, TrainAuc = mtrainauc, CalAUC = mcalauc)
  Multimodelauc <<- unique(rbind(Multimodelauc,mresultauc))
}

YTlogicregressionEVENCC <- glm(paste(target, ' ~ ', paste(select_variable_EVENcorelation, collapse=' + '), sep=''), data=YTtrain, family = binomial)

YTevalmodelauc(namea = 'logic regression for even cc',type = 'multivariate',name = 'logic regression', predmodelfortrain = 
                 predict(YTlogicregressionEVENCC, newdata=YTtrain, type = 'response'), truefortrain = YTtrain[, target],predmodelforcal = 
                 predict(YTlogicregressionEVENCC, newdata=YTCal, type = 'response'), trueforcal = YTCal[, target], plotname = "conma")

## Precision: 0.3510638 
## Recall: 0.6470588 
## Confusion matrix for logic regression for even cc 
## 
## 
## |      | [0,1.77e+07)| [1.77e+07,Inf]|
## |:-----|------------:|--------------:|
## |FALSE |           67|             61|
## |TRUE  |           18|             33|

YTevalmodelauc(namea = 'logic regression for even cc',type = 'multivariate',name = 'logic regression', predmodelfortrain = 
                 predict(YTlogicregressionEVENCC, newdata=YTtrain, type = 'response'), truefortrain = YTtrain[, target],predmodelforcal = 
                 predict(YTlogicregressionEVENCC, newdata=YTCal, type = 'response'), trueforcal = YTCal[, target], plotname = "loglike")

## Model Log Likelihood is -125.8907 and null Log Likelihood is -123.847 .

YTevalmodelauc(namea = 'logic regression for even cc',type = 'multivariate',name = 'logic regression', predmodelfortrain = 
                 predict(YTlogicregressionEVENCC, newdata=YTtrain, type = 'response'), truefortrain = YTtrain[, target],predmodelforcal = 
                 predict(YTlogicregressionEVENCC, newdata=YTCal, type = 'response'), trueforcal = YTCal[, target], plotname = "roc")

YTevalmodelauc(namea = 'logic regression for even cc',type = 'multivariate',name = 'logic regression', predmodelfortrain = 
                 predict(YTlogicregressionEVENCC, newdata=YTtrain, type = 'response'), truefortrain = YTtrain[, target],predmodelforcal = 
                 predict(YTlogicregressionEVENCC, newdata=YTCal, type = 'response'), trueforcal = YTCal[, target], plotname = "Predpplot")

YTevalmodelauc(namea = 'logic regression for even cc',type = 'multivariate',name = 'logic regression', predmodelfortrain = 
                 predict(YTlogicregressionEVENCC, newdata=YTtrain, type = 'response'), truefortrain = YTtrain[, target],predmodelforcal = 
                 predict(YTlogicregressionEVENCC, newdata=YTCal, type = 'response'), trueforcal = YTCal[, target], plotname = "Prereplot")

YTlogicregressionlog <- glm(paste(target, ' ~ ', paste(select_variable_log, collapse=' + '), sep=''), data=YTtrain, family = binomial)

YTevalmodelauc(namea = 'logic regression for log',type = 'multivariate',name = 'logic regression', predmodelfortrain = predict(YTlogicregressionlog, newdata=YTtrain, type = 'response'), truefortrain = YTtrain[, target],predmodelforcal = predict(YTlogicregressionlog, newdata=YTCal, type = 'response'), trueforcal = YTCal[, target], plotname = "conma")

## Precision: 0.9893617 
## Recall: 1 
## Confusion matrix for logic regression for log 
## 
## 
## |      | [0,1.77e+07)| [1.77e+07,Inf]|
## |:-----|------------:|--------------:|
## |FALSE |           85|              1|
## |TRUE  |            0|             93|

YTevalmodelauc(namea = 'logic regression for log',type = 'multivariate',name = 'logic regression', predmodelfortrain = predict(YTlogicregressionlog, newdata=YTtrain, type = 'response'), truefortrain = YTtrain[, target],predmodelforcal = predict(YTlogicregressionlog, newdata=YTCal, type = 'response'), trueforcal = YTCal[, target], plotname = "loglike")

## Model Log Likelihood is -36.83781 and null Log Likelihood is -123.847 .

YTevalmodelauc(namea = 'logic regression for log',type = 'multivariate',name = 'logic regression', predmodelfortrain = predict(YTlogicregressionlog, newdata=YTtrain, type = 'response'), truefortrain = YTtrain[, target],predmodelforcal = predict(YTlogicregressionlog, newdata=YTCal, type = 'response'), trueforcal = YTCal[, target], plotname = "roc")

YTevalmodelauc(namea = 'logic regression for log',type = 'multivariate',name = 'logic regression', predmodelfortrain = predict(YTlogicregressionlog, newdata=YTtrain, type = 'response'), truefortrain = YTtrain[, target],predmodelforcal = predict(YTlogicregressionlog, newdata=YTCal, type = 'response'), trueforcal = YTCal[, target], plotname = "Predpplot")

YTevalmodelauc(namea = 'logic regression for log',type = 'multivariate',name = 'logic regression', predmodelfortrain = predict(YTlogicregressionlog, newdata=YTtrain, type = 'response'), truefortrain = YTtrain[, target],predmodelforcal = predict(YTlogicregressionlog, newdata=YTCal, type = 'response'), trueforcal = YTCal[, target], plotname = "Prereplot")

YTlogicregressionODDCC <- glm(paste(target, ' ~ ', paste(select_variable_ODDcorelation, collapse=' + '), sep=''), data=YTtrain, family = binomial)

YTevalmodelauc(namea = 'logic regression for Odd cc',type = 'multivariate',name = 'logic regression', predmodelfortrain = predict(YTlogicregressionODDCC, newdata=YTtrain, type = 'response'), truefortrain = YTtrain[, target],predmodelforcal = predict(YTlogicregressionODDCC, newdata=YTCal, type = 'response'), trueforcal = YTCal[, target], plotname = "conma")

## Precision: 0.5744681 
## Recall: 0.7714286 
## Confusion matrix for logic regression for Odd cc 
## 
## 
## |      | [0,1.77e+07)| [1.77e+07,Inf]|
## |:-----|------------:|--------------:|
## |FALSE |           69|             40|
## |TRUE  |           16|             54|

YTevalmodelauc(namea = 'logic regression for Odd cc',type = 'multivariate',name = 'logic regression', predmodelfortrain = predict(YTlogicregressionODDCC, newdata=YTtrain, type = 'response'), truefortrain = YTtrain[, target],predmodelforcal = predict(YTlogicregressionODDCC, newdata=YTCal, type = 'response'), trueforcal = YTCal[, target], plotname = "loglike")

## Model Log Likelihood is -108.2383 and null Log Likelihood is -123.847 .

YTevalmodelauc(namea = 'logic regression for Odd cc',type = 'multivariate',name = 'logic regression', predmodelfortrain = predict(YTlogicregressionODDCC, newdata=YTtrain, type = 'response'), truefortrain = YTtrain[, target],predmodelforcal = predict(YTlogicregressionODDCC, newdata=YTCal, type = 'response'), trueforcal = YTCal[, target], plotname = "roc")

YTevalmodelauc(namea = 'logic regression for Odd cc',type = 'multivariate',name = 'logic regression', predmodelfortrain = predict(YTlogicregressionODDCC, newdata=YTtrain, type = 'response'), truefortrain = YTtrain[, target],predmodelforcal = predict(YTlogicregressionODDCC, newdata=YTCal, type = 'response'), trueforcal = YTCal[, target], plotname = "Predpplot")

YTevalmodelauc(namea = 'logic regression for Odd cc',type = 'multivariate',name = 'logic regression', predmodelfortrain = predict(YTlogicregressionODDCC, newdata=YTtrain, type = 'response'), truefortrain = YTtrain[, target],predmodelforcal = predict(YTlogicregressionODDCC, newdata=YTCal, type = 'response'), trueforcal = YTCal[, target], plotname = "Prereplot")

YTEVENCCDT <- rpart(paste(target, ' ~ ', paste(select_variable_EVENcorelation, collapse=' + '), sep=''),data=YTtrain)
rpart.plot(YTEVENCCDT)

YTevalmodelauc(namea = 'Decision Tree for even cc', type = 'multivariate', name = 'Decision Tree', predmodelfortrain = predict(YTEVENCCDT, newdata=YTtrain)[,pos.label], truefortrain = YTtrain[,target], predmodelforcal = predict(YTEVENCCDT, newdata=YTCal)[,pos.label],trueforcal = YTCal[,target], plotname = "conma")

## Precision: 0.5744681 
## Recall: 0.6352941 
## Confusion matrix for Decision Tree for even cc 
## 
## 
## |      | [0,1.77e+07)| [1.77e+07,Inf]|
## |:-----|------------:|--------------:|
## |FALSE |           54|             40|
## |TRUE  |           31|             54|

YTevalmodelauc(namea = 'Decision Tree for even cc', type = 'multivariate', name = 'Decision Tree', predmodelfortrain = predict(YTEVENCCDT, newdata=YTtrain)[,pos.label], truefortrain = YTtrain[,target], predmodelforcal = predict(YTEVENCCDT, newdata=YTCal)[,pos.label],trueforcal = YTCal[,target], plotname = "loglike")

## Model Log Likelihood is -123.9011 and null Log Likelihood is -123.847 .

YTevalmodelauc(namea = 'Decision Tree for even cc', type = 'multivariate', name = 'Decision Tree', predmodelfortrain = predict(YTEVENCCDT, newdata=YTtrain)[,pos.label], truefortrain = YTtrain[,target], predmodelforcal = predict(YTEVENCCDT, newdata=YTCal)[,pos.label],trueforcal = YTCal[,target], plotname = "roc")

YTevalmodelauc(namea = 'Decision Tree for even cc', type = 'multivariate', name = 'Decision Tree', predmodelfortrain = predict(YTEVENCCDT, newdata=YTtrain)[,pos.label], truefortrain = YTtrain[,target], predmodelforcal = predict(YTEVENCCDT, newdata=YTCal)[,pos.label],trueforcal = YTCal[,target], plotname = "Predpplot")

YTevalmodelauc(namea = 'Decision Tree for even cc', type = 'multivariate', name = 'Decision Tree', predmodelfortrain = predict(YTEVENCCDT, newdata=YTtrain)[,pos.label], truefortrain = YTtrain[,target], predmodelforcal = predict(YTEVENCCDT, newdata=YTCal)[,pos.label],trueforcal = YTCal[,target], plotname = "Prereplot")

YTlogDT <- rpart(paste(target, ' ~ ', paste(select_variable_log, collapse=' + '), sep=''),data=YTtrain)
rpart.plot(YTlogDT)

YTevalmodelauc(namea = 'Decision Tree for log', type = 'multivariate', name = 'Decision Tree', predmodelfortrain = predict(YTlogDT, newdata=YTtrain)[,pos.label], truefortrain = YTtrain[,target], predmodelforcal = predict(YTlogDT, newdata=YTCal)[,pos.label],trueforcal = YTCal[,target], plotname = "conma")

## Precision: 0.9893617 
## Recall: 1 
## Confusion matrix for Decision Tree for log 
## 
## 
## |      | [0,1.77e+07)| [1.77e+07,Inf]|
## |:-----|------------:|--------------:|
## |FALSE |           85|              1|
## |TRUE  |            0|             93|

YTevalmodelauc(namea = 'Decision Tree for log', type = 'multivariate', name = 'Decision Tree', predmodelfortrain = predict(YTlogDT, newdata=YTtrain)[,pos.label], truefortrain = YTtrain[,target], predmodelforcal = predict(YTlogDT, newdata=YTCal)[,pos.label],trueforcal = YTCal[,target], plotname = "loglike")

## Model Log Likelihood is -5.462418 and null Log Likelihood is -123.847 .

YTevalmodelauc(namea = 'Decision Tree for log', type = 'multivariate', name = 'Decision Tree', predmodelfortrain = predict(YTlogDT, newdata=YTtrain)[,pos.label], truefortrain = YTtrain[,target], predmodelforcal = predict(YTlogDT, newdata=YTCal)[,pos.label],trueforcal = YTCal[,target], plotname = "roc")

YTevalmodelauc(namea = 'Decision Tree for log', type = 'multivariate', name = 'Decision Tree', predmodelfortrain = predict(YTlogDT, newdata=YTtrain)[,pos.label], truefortrain = YTtrain[,target], predmodelforcal = predict(YTlogDT, newdata=YTCal)[,pos.label],trueforcal = YTCal[,target], plotname = "Predpplot")

YTevalmodelauc(namea = 'Decision Tree for log', type = 'multivariate', name = 'Decision Tree', predmodelfortrain = predict(YTlogDT, newdata=YTtrain)[,pos.label], truefortrain = YTtrain[,target], predmodelforcal = predict(YTlogDT, newdata=YTCal)[,pos.label],trueforcal = YTCal[,target], plotname = "Prereplot")

YTODDCCDT <- rpart(paste(target, ' ~ ', paste(select_variable_ODDcorelation, collapse=' + '), sep=''),data=YTtrain)
rpart.plot(YTODDCCDT)

YTevalmodelauc(namea = 'Decision Tree for odd cc', type = 'multivariate', name = 'Decision Tree', predmodelfortrain = predict(YTODDCCDT, newdata=YTtrain)[,pos.label], truefortrain = YTtrain[,target], predmodelforcal = predict(YTODDCCDT, newdata=YTCal)[,pos.label],trueforcal = YTCal[,target], plotname = "conma")

## Precision: 0.5319149 
## Recall: 0.7142857 
## Confusion matrix for Decision Tree for odd cc 
## 
## 
## |      | [0,1.77e+07)| [1.77e+07,Inf]|
## |:-----|------------:|--------------:|
## |FALSE |           65|             44|
## |TRUE  |           20|             50|

YTevalmodelauc(namea = 'Decision Tree for odd cc', type = 'multivariate', name = 'Decision Tree', predmodelfortrain = predict(YTODDCCDT, newdata=YTtrain)[,pos.label], truefortrain = YTtrain[,target], predmodelforcal = predict(YTODDCCDT, newdata=YTCal)[,pos.label],trueforcal = YTCal[,target], plotname = "loglike")

## Model Log Likelihood is -114.8999 and null Log Likelihood is -123.847 .

YTevalmodelauc(namea = 'Decision Tree for odd cc', type = 'multivariate', name = 'Decision Tree', predmodelfortrain = predict(YTODDCCDT, newdata=YTtrain)[,pos.label], truefortrain = YTtrain[,target], predmodelforcal = predict(YTODDCCDT, newdata=YTCal)[,pos.label],trueforcal = YTCal[,target], plotname = "roc")

YTevalmodelauc(namea = 'Decision Tree for odd cc', type = 'multivariate', name = 'Decision Tree', predmodelfortrain = predict(YTODDCCDT, newdata=YTtrain)[,pos.label], truefortrain = YTtrain[,target], predmodelforcal = predict(YTODDCCDT, newdata=YTCal)[,pos.label],trueforcal = YTCal[,target], plotname = "Predpplot")

YTevalmodelauc(namea = 'Decision Tree for odd cc', type = 'multivariate', name = 'Decision Tree', predmodelfortrain = predict(YTODDCCDT, newdata=YTtrain)[,pos.label], truefortrain = YTtrain[,target], predmodelforcal = predict(YTODDCCDT, newdata=YTCal)[,pos.label],trueforcal = YTCal[,target], plotname = "Prereplot")

YTknntrainevenCC <- knn(YTtrain[select_variable_EVENcorelation], YTtrain[select_variable_EVENcorelation], YTtrain[,target], k=3, prob=T)
YTknnCalevenCC <- knn(YTtrain[select_variable_EVENcorelation], YTCal[select_variable_EVENcorelation], YTtrain[,target], k=3, prob=T)
YTtrainknnProbevenCC <- attributes(YTknntrainevenCC)$prob
YTcalknnProbevenCC <- attributes(YTknnCalevenCC)$prob

knntable <- table(actual=YTCal[,target], predicted=YTknnCalevenCC)
      
TP <- knntable[2, 2]
FP <- knntable[1, 2]
FN <- knntable[2, 1]
      
precision <- TP / (TP + FP)
recall <- TP / (TP + FN)
      
cat("Precision:", precision, "\nRecall:", recall, "\nConfusion matrix for knn mode.\n")

## Precision: 0.5494505 
## Recall: 0.5319149 
## Confusion matrix for knn mode.

print(kable(knntable))

## 
## 
## |               | [0,1.77e+07)| [1.77e+07,Inf]|
## |:--------------|------------:|--------------:|
## |[0,1.77e+07)   |           44|             41|
## |[1.77e+07,Inf] |           44|             50|

YTKNNcalauc <- YTAUC(YTcalknnProbevenCC, YTCal[,target])
YTknntrainevenCCauc <- YTAUC(YTtrainknnProbevenCC, YTtrain[,target])
KNNauc <- data.frame(Pred = "KNN for even CC", Type = 'multivariate', TrainAuc = YTknntrainevenCCauc, CalAUC = YTKNNcalauc)
Multimodelauc <- unique(rbind(Multimodelauc,KNNauc))
kable(KNNauc)

knnplot <- function(predmodel, trueV, titleString="ROC plot") {
      perf <- performance(prediction(predmodel, trueV), "tpr", "fpr")
      YTROknn <- data.frame(
        FPR = unlist(perf@x.values),
        TPR = unlist(perf@y.values),
        threshold = unlist(perf@alpha.values)
      )
      
      ggplot(YTROknn, aes(x = FPR, y = TPR)) +
        geom_line(color = "blue") +
        geom_abline(intercept = 0, slope = 1, linetype = 'dashed', color = "red") +
        labs(x = "FPR", y = "TPR", title=titleString) +
        theme_minimal()
}
knnplot(YTcalknnProbevenCC , YTCal[,target], titleString="AUC of kNN predictions for even CC")

knnprreplot <- function(predmodel, trueV, namea) {          
  YTprec <- performance(prediction(predmodel, trueV), measure="prec")
          YTrec <- performance(prediction(predmodel, trueV), measure="rec")
          precision <- (YTprec@y.values)[[1]]
          YTxprec <- (YTprec@x.values)[[1]] 
          recall <- (YTrec@y.values)[[1]]
          rocFrame <- data.frame(threshold=YTxprec, precision=precision, recall=recall)
          ggplot(rocFrame, aes(x=threshold)) +
              geom_line(aes(y=precision, color="precision")) +
              geom_line(aes(y=recall, color="recall")) +
              coord_cartesian(xlim = c(0, 1), ylim = c(0, 1)) +
              labs(x = "threshold", y = "percentage", title = paste("precision and recall plot for", namea)) + 
              theme_minimal()
          
}
knnprreplot(YTcalknnProbevenCC , YTCal[,target], 'KNN for even cc')

YTknntrainlog <- knn(YTtrain[select_variable_log], YTtrain[select_variable_log], YTtrain[,target], k=3, prob=T)
YTknnCal <- knn(YTtrain[select_variable_log], YTCal[select_variable_log], YTtrain[,target], k=3, prob=T)
YTtrainknnProblog <- attributes(YTknntrainlog)$prob
YTcalknnProblog <- attributes(YTknnCal)$prob

knntable <- table(actual=YTCal[,target], predicted=YTknnCal)
      
TP <- knntable[2, 2]
FP <- knntable[1, 2]
FN <- knntable[2, 1]
      
precision <- TP / (TP + FP)
recall <- TP / (TP + FN)
      
cat("Precision:", precision, "\nRecall:", recall, "\nConfusion matrix for knn mode.\n")

## Precision: 1 
## Recall: 0.9893617 
## Confusion matrix for knn mode.

print(kable(knntable))

## 
## 
## |               | [0,1.77e+07)| [1.77e+07,Inf]|
## |:--------------|------------:|--------------:|
## |[0,1.77e+07)   |           85|              0|
## |[1.77e+07,Inf] |            1|             93|

YTKNNcalauc <- YTAUC(YTcalknnProblog, YTCal[,target])
YTKNNtrainauc <- YTAUC(YTtrainknnProblog, YTtrain[,target])
KNNauc <- data.frame(Pred = "KNN", Type = 'multivariate', TrainAuc = YTKNNtrainauc, CalAUC = YTKNNcalauc)
Multimodelauc <- unique(rbind(Multimodelauc,KNNauc))
kable(KNNauc)

knnplot <- function(predmodel, trueV, titleString="ROC plot for log") {
      perf <- performance(prediction(predmodel, trueV), "tpr", "fpr")
      YTROknn <- data.frame(
        FPR = unlist(perf@x.values),
        TPR = unlist(perf@y.values),
        threshold = unlist(perf@alpha.values)
      )
      
      ggplot(YTROknn, aes(x = FPR, y = TPR)) +
        geom_line(color = "blue") +
        geom_abline(intercept = 0, slope = 1, linetype = 'dashed', color = "red") +
        labs(x = "FPR", y = "TPR", title=titleString) +
        theme_minimal()
}
knnplot(YTcalknnProblog , YTCal[,target], titleString="AUC of kNN predictions for log")

knnprreplot <- function(predmodel, trueV, namea) {          
  YTprec <- performance(prediction(predmodel, trueV), measure="prec")
          YTrec <- performance(prediction(predmodel, trueV), measure="rec")
          precision <- (YTprec@y.values)[[1]]
          YTxprec <- (YTprec@x.values)[[1]] 
          recall <- (YTrec@y.values)[[1]]
          rocFrame <- data.frame(threshold=YTxprec, precision=precision, recall=recall)
          ggplot(rocFrame, aes(x=threshold)) +
              geom_line(aes(y=precision, color="precision")) +
              geom_line(aes(y=recall, color="recall")) +
              coord_cartesian(xlim = c(0, 1), ylim = c(0, 1)) +
              labs(x = "threshold", y = "percentage", title = paste("precision and recall plot for", namea)) + 
              theme_minimal()
          
}
knnprreplot(YTcalknnProblog , YTCal[,target], 'KNN for log')

YTknntrainODDcc <- knn(YTtrain[select_variable_ODDcorelation], YTtrain[select_variable_ODDcorelation], YTtrain[,target], k=3, prob=T)
YTknnCalODdcc <- knn(YTtrain[select_variable_ODDcorelation], YTCal[select_variable_ODDcorelation], YTtrain[,target], k=3, prob=T)
YTtrainknnProbodCC <- attributes(YTknntrainODDcc)$prob
YTcalknnProbodCC <- attributes(YTknnCalODdcc)$prob

knntable <- table(actual=YTCal[,target], predicted=YTknnCalODdcc)
      
TP <- knntable[2, 2]
FP <- knntable[1, 2]
FN <- knntable[2, 1]
      
precision <- TP / (TP + FP)
recall <- TP / (TP + FN)
      
cat("Precision:", precision, "\nRecall:", recall, "\nConfusion matrix for knn mode.\n")

## Precision: 0.6292135 
## Recall: 0.5957447 
## Confusion matrix for knn mode.

print(kable(knntable))

## 
## 
## |               | [0,1.77e+07)| [1.77e+07,Inf]|
## |:--------------|------------:|--------------:|
## |[0,1.77e+07)   |           52|             33|
## |[1.77e+07,Inf] |           38|             56|

YTKNNcalauc <- YTAUC(YTcalknnProbodCC, YTCal[,target])
YTKNNtrainauc <- YTAUC(YTtrainknnProbodCC, YTtrain[,target])
KNNauc <- data.frame(Pred = "KNN", Type = 'multivariate', TrainAuc = YTKNNtrainauc, CalAUC = YTKNNcalauc)
Multimodelauc <- unique(rbind(Multimodelauc,KNNauc))
kable(KNNauc)

knnplot <- function(predmodel, trueV, titleString="ROC plot") {
      perf <- performance(prediction(predmodel, trueV), "tpr", "fpr")
      YTROknn <- data.frame(
        FPR = unlist(perf@x.values),
        TPR = unlist(perf@y.values),
        threshold = unlist(perf@alpha.values)
      )
      
      ggplot(YTROknn, aes(x = FPR, y = TPR)) +
        geom_line(color = "blue") +
        geom_abline(intercept = 0, slope = 1, linetype = 'dashed', color = "red") +
        labs(x = "FPR", y = "TPR", title=titleString) +
        theme_minimal()
}
knnplot(YTcalknnProbodCC , YTCal[,target], titleString="AUC of kNN predictions of odd cc")

knnprreplot <- function(predmodel, trueV, namea) {          
  YTprec <- performance(prediction(predmodel, trueV), measure="prec")
          YTrec <- performance(prediction(predmodel, trueV), measure="rec")
          precision <- (YTprec@y.values)[[1]]
          YTxprec <- (YTprec@x.values)[[1]] 
          recall <- (YTrec@y.values)[[1]]
          rocFrame <- data.frame(threshold=YTxprec, precision=precision, recall=recall)
          ggplot(rocFrame, aes(x=threshold)) +
              geom_line(aes(y=precision, color="precision")) +
              geom_line(aes(y=recall, color="recall")) +
              coord_cartesian(xlim = c(0, 1), ylim = c(0, 1)) +
              labs(x = "threshold", y = "percentage", title = paste("precision and recall plot for", namea)) + 
              theme_minimal()
          
}
knnprreplot(YTcalknnProbodCC , YTCal[,target], 'KNN of odd cc')

kable(Multimodellikelog)

colnames(AUCYTNV) <- colnames(Multimodelauc)
AUCtotaltable <- rbind(Multimodelauc, AUCYTNV)
kable(head(AUCtotaltable[order(AUCtotaltable$CalAUC, decreasing = T), ], 10))

YTdiffModel <- list(
  list(predmodel = predict(YTlogicregressionEVENCC, newdata = YTtrain, type = 'response'), trueV = YTtrain[, target], namea = "LogisticRegression for evencc"),
  list(predmodel = predict(YTlogicregressionODDCC, newdata = YTtrain, type = 'response'), trueV = YTtrain[, target], namea = "LogisticRegression for ODDcc"),
  list(predmodel = predict(YTlogicregressionEVENCC, newdata = YTtrain, type = 'response'), trueV = YTtrain[, target], namea = "LogisticRegression for evencc"),
  list(predmodel = predict(YTEVENCCDT, newdata = YTtrain)[, pos.label], trueV = YTtrain[, target], namea = "DecisionTree for even CC"),
  list(predmodel = predict(YTlogDT, newdata = YTtrain)[, pos.label], trueV = YTtrain[, target], namea = "DecisionTree for log"),
  list(predmodel = predict(YTODDCCDT, newdata = YTtrain)[, pos.label], trueV = YTtrain[, target], namea = "DecisionTree for odd CC"),
  list(predmodel = YTcalknnProbevenCC, trueV = YTCal[, target], namea = "KNN for even cc"),
  list(predmodel = YTcalknnProblog, trueV = YTCal[, target], namea = "KNN for log"),
  list(predmodel = YTcalknnProbodCC, trueV = YTCal[, target], namea = "KNN for odd cc")
)

YTmodels <- do.call(rbind, lapply(YTdiffModel, as.data.frame))

YTrocG <- list()

for (model_name in unique(YTmodels$namea)) {
  model_subset <- YTmodels[YTmodels$namea == model_name, ]
  perf <- performance(prediction(model_subset$predmodel, model_subset$trueV), "tpr", "fpr")
  YTROCMulti <- data.frame(
    FPR = unlist(perf@x.values),
    TPR = unlist(perf@y.values),
    threshold = unlist(perf@alpha.values),
    namea = model_name
  )
  YTrocG[[model_name]] <- YTROCMulti
}

combined_roc <- do.call(rbind, YTrocG)

ggplot(combined_roc, aes(x = FPR, y = TPR, color = namea)) +
  geom_line() +
  geom_abline(intercept = 0, slope = 1, linetype = 'dashed', color = "red") + 
  labs(x = "FPR", y = "TPR", title = "ROC Curves for Multiple Models") +
  labs(color = "Model") +
  theme_minimal()

Final select model

As a result, Logistic Regression of Correlation Coefficient for Odd Variables is selected as the final mode. The evaluation of the final selected model by usign the testing set is shown below.

predicted_values <- predict(YTlogicregressionODDCC, newdata=YTtest, type = 'response')
YTmodeleval('logistic regression fpr for final model', predicted_values, YTtest[, target], do.print = TRUE, selection = "conma")

## Precision: 0.5841584 
## Recall: 0.8082192 
## Confusion matrix for logistic regression fpr for final model 
## 
## 
## |      | [0,1.77e+07)| [1.77e+07,Inf]|
## |:-----|------------:|--------------:|
## |FALSE |           88|             42|
## |TRUE  |           14|             59|

## [1] 0.7639293

predicted_values <- predict(YTlogicregressionODDCC, newdata=YTtest, type = 'response')
YTmodeleval('logistic regression fpr for final model', predicted_values, YTtest[, target], do.print = TRUE, selection = "loglike")

## Model Log Likelihood is -115.3738 and null Log Likelihood is -140.7064 .

## [1] 0.7639293

predicted_values <- predict(YTlogicregressionODDCC, newdata=YTtest, type = 'response')
YTmodeleval('logistic regression fpr for final model', predicted_values, YTtest[, target], do.print = TRUE, selection = "roc")

## [1] 0.7639293

predicted_values <- predict(YTlogicregressionODDCC, newdata=YTtest, type = 'response')
YTmodeleval('logistic regression fpr for final model', predicted_values, YTtest[, target], do.print = TRUE, selection = "Predpplot")

## [1] 0.7639293

predicted_values <- predict(YTlogicregressionODDCC, newdata=YTtest, type = 'response')
YTmodeleval('logistic regression fpr for final model', predicted_values, YTtest[, target], do.print = TRUE, selection = "Prereplot")

## [1] 0.7639293

Conclusion

In terms of model selection, single variable model, multivariate models, logistic regression, decision tree, and KNN model are used. Through the comparison of AUC, confusion matrix, ROC curve, predicted distribution, and precision and recall plot, logistic regression is the most accurate model among the other model. Besides, log likelihood, odd variables, even variables are used and odd variables provides more accurate predictions. Hence, logistic regression with correlation coefficient of odd variables is chosen as the final model to do the predictions.

Data Preparation

YT_cluster <- read.csv("./youtube_UTF_8.csv", header = TRUE, encoding = "UTF-8")
YTNV_cluster <- YT_cluster %>%
  mutate(
    allunemployment = (Unemployment.rate / 100) * Population,
    unemployment_urban = (Unemployment.rate / 100) * Urban_population,
    educationall = (Gross.tertiary.education.enrollment.... / 100) * Population,
    educationurban = (Gross.tertiary.education.enrollment.... / 100) * Urban_population,
    averagemonthearning = (lowest_yearly_earnings+ highest_yearly_earnings) / 2,
    averageyearearning = (lowest_monthly_earnings+ highest_monthly_earnings) / 2
  ) %>%
  select(-highest_monthly_earnings, -lowest_monthly_earnings, -lowest_yearly_earnings, -highest_yearly_earnings, -Gross.tertiary.education.enrollment...., -Unemployment.rate, -video_views_rank, -country_rank, -channel_type_rank, -created_year, -created_date) %>%
  filter(rowSums(is.na(select(., (ncol(.)-9):(ncol(.)-2)))) == 0) %>%
  select_if(is.numeric)  %>%
  mutate(across(everything(), ~ ifelse(is.na(.), mean(., na.rm = TRUE), .)))

Descriptive statistics for the new data like minimum, median, mean, standard deviation, and maximum are performed for data checking.

YTNV_cluster_summary <- data.frame(
  Min = sapply(YTNV_cluster, min),
  Med = sapply(YTNV_cluster, median),
  Mean = sapply(YTNV_cluster, mean),
  Max = sapply(YTNV_cluster, max),
  SD = sapply(YTNV_cluster, sd)
)
kable(YTNV_cluster_summary, align = "lccccc")

Data Scaling

Standardize the data to have a mean of 0 and a standard deviation of 1 for each variable. This makes sure that no single variable’s scale improperly affects the cluster formation.

vars.to.use <- colnames(YTNV_cluster)[-2]
scaled_df <- scale(YTNV_cluster[,vars.to.use])
attr(scaled_df, "scaled:center")

##                             rank                      video.views 
##                     5.052385e+02                     1.123485e+10 
##                          uploads video_views_for_the_last_30_days 
##                     1.040791e+04                     1.823042e+08 
##     subscribers_for_last_30_days                       Population 
##                     3.519155e+05                     4.303873e+08 
##                 Urban_population                         Latitude 
##                     2.242150e+08                     2.663278e+01 
##                        Longitude                  allunemployment 
##                    -1.412815e+01                     3.491125e+07 
##               unemployment_urban                     educationall 
##                     2.187965e+07                     2.031631e+08 
##                   educationurban              averagemonthearning 
##                     1.305702e+08                     4.002513e+06 
##               averageyearearning 
##                     3.333274e+05

attr(scaled_df, "scaled:scale")

##                             rank                      video.views 
##                     2.873427e+02                     1.479528e+10 
##                          uploads video_views_for_the_last_30_days 
##                     3.631451e+04                     4.220909e+08 
##     subscribers_for_last_30_days                       Population 
##                     5.061499e+05                     4.727947e+08 
##                 Urban_population                         Latitude 
##                     1.546874e+08                     2.056053e+01 
##                        Longitude                  allunemployment 
##                     8.476081e+01                     2.670902e+07 
##               unemployment_urban                     educationall 
##                     1.593357e+07                     1.403330e+08 
##                   educationurban              averagemonthearning 
##                     9.038184e+07                     7.439542e+06 
##               averageyearearning 
##                     6.194923e+05

Hierarchical Clustering

hc <- hclust(dist(scale(YTNV_cluster[,colnames(YTNV_cluster)[-2]]), method = "euclidean"), method = "ward.D2")

hcplot <- as.dendrogram(hc)
plot(hcplot, ylab = "Height", leaflab = "none", horiz = FALSE)
hightnu <- c(61.5, 55.4, 33.7)
coline <- c("royalblue2", "royalblue3", "royalblue4")  
for (i in 1:length(hightnu)) {
  abline(h = hightnu[i], col = coline[i], lty = 8, lwd = 1)
}
text(3, 20, paste("k >=", 5), col = "blue", cex = 0.6)
text(3, c(59, 48), paste("k =", 3:4), col = "blue", cex = 0.6)
text(3, 65, paste("k =", 2), col = "blue", cex = 0.6)

A dendogram illustrates the structure of a dataset by displaying the the link between clusters of data points. A new cluster is created when two branches combined. The blue dotted horizontal lines in the dendogram reflect various levels at which clusters form, and k-values show overall there are how many clusters. Based on the dendogram above, two major clusters are formed when k is equal to 2. The three and four clusters are further divided at k=3 and k=4. At k greater than and equal to 5, the dendogram predicts further subdivisions, but the dendogram above does not explicitly display them. Hence, we need to find the suitable k value.

wss <- function(cluster) {
  c0 <- colMeans(cluster)
  sum(apply(cluster, 1, function(row) sum((row - c0) ^ 2)))
}

# Calculate the total WSS
Totalwss <- function(df, labels) {
  sumwss <- 0
  for (i in 1:length(unique(labels))) {
    sumwss <- sumwss + wss(subset(df, labels == i))
  }
  sumwss
}

indexch <- function(df, kmax, data) {
  npts <- nrow(df)
  wss.value <- numeric(kmax)
  wss.value[1] <- wss(df)
  for (k in 2:kmax) {
    labels <- cutree(data, k)
    wss.value[k] <- Totalwss(df, labels)
  }
  B <- (wss(df) - wss.value)/ (0:(kmax - 1))
  w <- wss.value / (npts - 1:kmax)
  data.frame(k = 1:kmax, CH_index = B/w, WSS = wss.value)
}

The CH index is a benchmark for figuring out the ideal cluster count in hierarchical clustering. It strikes a balance between two factors of clustering: Compactness (how far apart elements in a cluster are from one another) and Separation. By determining the distance between data points within a cluster and its center, compactness is determined by how tightly a cluster is. A smaller number indicates clearly defined cluster. On the other hand, separation measures the separation between clusters by computing the distance between the cluster centers and the centroid of the entire dataset. A higher value indicates clear, separated clusters. Hence, a good clustering should have a higher CH index that signifies well-separated and have little internal variation. The k value with the highest CH index that denotes the most distinct and consistent clustering structure, is the goal to identify the ideal cluster number.

ch_criterion <- indexch(scaled_df, 10, hc)
grid.arrange(
  ggplot(ch_criterion, aes(x = k, y = CH_index)) +
    geom_point() +
    geom_line(colour = "red") +
    scale_x_continuous(breaks = 1:10, labels = 1:10) +
    labs(y = "CH index") +
    theme_minimal(),
  ggplot(ch_criterion, aes(x = k, y = WSS), color = "blue") +
    geom_point() + geom_line(colour = "blue") +
    scale_x_continuous(breaks = 1:10, labels = 1:10) +
    theme_minimal(),
  ncol = 2
)

Based on the plot, CH index is maximized at k=4 that appears to show a large peak.

plot(hcplot, ylab = "Height", leaflab = "none", horiz = F)
rect.hclust(hc, k = 4)

The dendogram shows the formation of four clusters by drawing red rectangles around the clusters.

kmean Clusters

kmch <- kmeansruns(scaled_df, krange=1:10, criterion="ch")
kmasw <- kmeansruns(scaled_df, krange=1:10, criterion="asw")

kmCritframe <- data.frame(k=1:10, ch=kmch$crit,asw=kmasw$crit)
grid.arrange(ggplot(kmCritframe, aes(x=k, y=ch)) +
  geom_point() + geom_line(colour="red") +
  scale_x_continuous(breaks=1:10, labels=1:10) +
  labs(y="CH index") +
    theme_minimal(),
ggplot(kmCritframe, aes(x=k, y=asw)) +
  geom_point() + geom_line(colour="blue") +
  scale_x_continuous(breaks=1:10, labels=1:10) +
  labs(y="ASW") +
    theme_minimal(),
nrow=1)

Plot above demonstrates peaks at particular k values, which may point to the best possible grouping. The best cluster count, according to the CH index plot, is k = 3, while the optimum cluster count, according to the ASW plot is k = 4.

Clustersing Visualising

princ <- prcomp(scaled_df)
project2D <- as.data.frame(predict(princ, newdata = scaled_df)[, 1:2])

convexHullcheck <- function(proj2ddata, groups) {
  do.call(rbind, lapply(unique(groups), function(c) {
    f <- subset(proj2ddata, cluster == c)
    f[chull(f), ]
  }))
}

plot_clusters <- function(number) {
  groups <- cutree(hc, number)
  convexHulldf <<- cbind(
    project2D,
    cluster = as.factor(groups),
    subscriber = YTNV_cluster$subscribers
  )
  convex_hull <<- convexHullcheck(convexHulldf, groups)
  ggplot(convexHulldf, aes(x = PC1, y = PC2)) +
    geom_point(aes(shape = cluster, color = cluster, alpha = 0.1)) +
    geom_polygon(data = convex_hull, aes(group = cluster, fill = cluster), alpha = 0.4, linetype = 0) +
    labs(title = sprintf("k = %d", number)) +
    theme(legend.position = "none") +
    theme_minimal()
}

# Generate plots for k=3 and k=4
plots <- lapply(c(3,4), plot_clusters)

# Display plots side by side
grid.arrange(plots[[1]], plots[[2]], ncol = 2)

The data is divided into two and three clusters and decide which cluster is more ideal.

fviz_cluster(kmeans(scaled_df, 2), scaled_df)

fviz_cluster(kmeans(scaled_df, 3), scaled_df)

It can be observed when k = 2 the clusters are clearly separated from one another, however when k = 3 there the cluster is overlapping.

Extract members of each cluster

At this point, the dataset will be clustered using the k-mean algorithm, and the criteria for Cluster 1 and Cluster 2 will be determined based on the traits of each cluster after clustering.

set.seed(4009)
kM <- kmeans(scaled_df, 2)
CLustersummary <- cbind(YTNV_cluster, cluster = kM$cluster)
clusters <- table(CLustersummary$cluster)

kable(table(clusters))

kmresult <- kmeans(scaled_df, 2, nstart=100, iter.max=100)
kmresult$centers

##          rank video.views       uploads video_views_for_the_last_30_days
## 1 -0.06342256   0.0784036  0.0006477713                       0.02517707
## 2  0.07838378  -0.0968988 -0.0008005789                      -0.03111628
##   subscribers_for_last_30_days Population Urban_population   Latitude
## 1                   0.01698760  0.5539946        0.7594241  0.2288399
## 2                  -0.02099493 -0.6846805       -0.9385703 -0.2828227
##    Longitude allunemployment unemployment_urban educationall educationurban
## 1 -0.2393459       0.8264962          0.8053987    0.8560584      0.7907666
## 2  0.2958070      -1.0214646         -0.9953902   -1.0580003     -0.9773064
##   averagemonthearning averageyearearning
## 1          0.05044994         0.05064612
## 2         -0.06235096        -0.06259341

kmresult$size

## [1] 482 390

Considering the prior analysis, the dataset is divided into 2 sets of 872 data. Cluster 1 has been assigned 482 data whereas Cluster 2 has been assigned 390 data. The data points in Cluster 1 appear to have relatively higher values for the variables video.views, uploads, urban_population, and education. On the other hand, Cluster 2 appears to reflect data points with average values that are lower for the majority of the variables.

HCdata <- clusterboot(scaled_df, clustermethod=hclustCBI,method="ward.D2", k=2)

summary(HCdata$result)

##               Length Class  Mode     
## result          7    hclust list     
## noise           1    -none- logical  
## nc              1    -none- numeric  
## clusterlist     2    -none- list     
## partition     872    -none- numeric  
## clustermethod   1    -none- character
## nccl            1    -none- numeric

1 - HCdata$bootbrd/100

## [1] 1 1

Both the Cluster 1 and Cluster 2 are very stable.

Conclusion

In a nutshell, a deeper comprehension of the underlying structures in the YouTube dataset has been made by clustering analysis. The ideal number of clusters is determined and k=4 clusters are the most effective for our dataset. Key ideas provide the foundation for the rationale behind this.

1. The testing indicates that four clusters provide the optimal clustering, which corresponds to the illustration above that uses the CH index to determine the ideal number pf clusters. This indicates that the data points are distinctly divided from the other groups and densely packed within their respective clusterings.

2. From the perspective of kmeans clustering, k=2 is preferred as it shows a clearer separation with no overlapping.

3. Finally, since the cluster plot in Kmean cluster is clearer and more dispersed than Hierarchical Clustering, and we can clearly figure out the clusters. We choose Kmean clustering.